Show the code
# Load libraries
library(tidyverse)
library(tidybayes)
library(bayesplot)
library(rstanarm)
library(here)
library(plotly)
library(sf)This document uses tidybayes to extract and visualize model fit params for the NEON data/ autotrophic uptake project
First created 2/25/2025 by Christa Torrens, with ongoing amendments
# Load libraries
library(tidyverse)
library(tidybayes)
library(bayesplot)
library(rstanarm)
library(here)
library(plotly)
library(sf)This keeps all of the core plotting code in 1 place, and allows site-level customization (mostly dataframes and titles)
##### Model fit - Npred + Nobs over time
plot_Nmodelfit_ts <- function(data, start_jday, end_jday, plot_title) {
plot <- data %>%
filter(model_jday >= start_jday & model_jday <= end_jday) %>%
ggplot(aes(x = mod_hours)) +
geom_point(aes(y = N_conc)) +
geom_line(aes(y = N_pred), col = 'red') +
labs(
x = "Time (h)",
y = expression("N" ~ (mmol ~ m^-3))
) +
ggtitle(plot_title) +
facet_wrap(~yr_jday) +
theme_bw()
}
plot_Nmodelfit_tsCI <- function(data, start_jday, end_jday, plot_title) {
plot <- data %>%
filter(model_jday >= start_jday & model_jday <= end_jday) %>%
ggplot(aes(x = mod_hours)) +
geom_ribbon(aes(ymin=N_pred.lower, ymax=N_pred.upper), alpha = 0.7, fill = "slategray") +
geom_point(aes(y = N_conc)) +
geom_line(aes(y = N_pred), col = 'red') +
labs(
x = "Time (h)",
y = expression("N" ~ (mmol ~ m^-3))
) +
ggtitle(plot_title) +
facet_wrap(~yr_jday) +
theme_bw()
}
##### Model fit - Npred v Nobs w. 1:1 line
plot_Npred_v_Nobs <- function(data, plot_title) {
plot <- ggplot(data=data, aes(x=N_conc, y=N_pred)) +
geom_point() +
labs( x=expression("measured N"~(mmol~m^-3)), # OR (mu*mol~L^-1)
y=expression("predicted N"~(mmol~m^-3)) # ~ = a space, * = no space
) +
ggtitle(plot_title) +
geom_abline(intercept = 0, slope = 1, color = "red") +
theme_bw()
}
##### Equifinality checks
# added error bars to all equifinality checks, 4/16/25
## K vs U
plot_KvU <- function(data, plot_title) {
plot <- ggplot(data=data, aes(y=K, x=U)) +
geom_errorbar(aes(ymin = K.lower, ymax = K.upper),
width = 0.2, color='grey70') + # CI bars
geom_point() +
labs(y=expression("modeled K"~(d^-1)),
x=expression("modeled U"~(mmol~m^-2~d^-1))
) +
ggtitle(plot_title) +
# geom_abline(intercept = 0, slope = 1, color = "red") + # doesn't show up - off the charts!
theme_bw()
}
## K vs N_e
plot_KvNe <- function(data, plot_title) {
plot <- ggplot(data=data, aes(y=K, x=N_e)) +
geom_errorbar(aes(ymin = K.lower, ymax = K.upper),
width = 0.2, color='grey70') + # CI bars
geom_point() +
labs(y=expression("modeled K"~(d^-1)),
x=expression("modeled N_e"~(mmol~m^-3))
) +
ggtitle(plot_title) +
# geom_abline(intercept = 0, slope = 1, color = "red") + # doesn't show up - off the charts!
theme_bw()
}
## U vs N_e
# Includes CI bars because this relationship tends to look a bit linear - checking whether that holds even w ci values in place
plot_UvNe <- function(data, plot_title) {
plot <- ggplot(data=data, aes(y=U, x=N_e)) +
geom_errorbar(aes(ymin = U.lower, ymax = U.upper),
width = 0.2, color='grey70') + # CI bars
geom_point() +
labs(y=expression("modeled U"~(mmol~m^-2~d^-1)),
x=expression("modeled N_e"~(mmol~m^-3))
) +
ggtitle(plot_title) +
theme_bw()
}
##### K over time
plot_K_ts <- function(data, plot_title) {
plot <- ggplot(data=data, aes(x=d, y=K)) +
geom_errorbar(aes(ymin = K.lower, ymax = K.upper),
width = 0.2, color='grey70') + # CI bars
geom_point() +
# geom_point(y = sumlight.real, color = 'gold') +
labs( x= "Index",
y=expression("modeled K"~(d^-1)) # ~ = a space, * = no space
) +
ggtitle(plot_title) +
theme_bw()
}
##### U over time
plot_U_ts <- function(data, plot_title) {
ggplot(data=data, aes(x=d, y=U)) +
geom_errorbar(aes(ymin = U.lower, ymax = U.upper),
width = 0.2, color='grey70') + # CI bars
geom_point() +
# geom_point(y = sumlight.real_wlou, color = 'gold') +
# ADD IN HIGH AND LOW CIs
labs( x= "Index",
y=expression("modeled U"~(mmol~m^-2~d^-1)) # ~ = a space, * = no space
) +
ggtitle(plot_title) +
theme_bw()
}
##### U vs. real sumlight
plot_UvLight <- function(data, plot_title) {
plot <- ggplot(data=data, aes(x=sumlightReal, y=U)) +
geom_point() +
labs( x=expression("measured daily light"~(w~m^-2)),
y=expression("modeled U"~(mmol~m^-2~d^-1)) # ~ = a space, * = no space
) +
scale_x_log10() + scale_y_log10() +
ggtitle(plot_title) +
theme_bw()
}
##### U + daily (or hourly?) sumlight over time
##### Autotrophic U vs total U
plot_autoUperc <- function(data, plot_title) {
plot <- ggplot(data=data, aes(x=d, y=U_auto_perc)) +
geom_point() +
xlab("Index") + ylab("Percent autotrophic uptake") +
ggtitle(plot_title) +
theme_bw()
}Outputs and visualizations for Upper Big Creek, Fresno, CA
########## Load model fit ______________________________________________________
bigc.fit <- readRDS(here("N_uptake_NEON/data/model_fits/bigc.fit.rds"))
# if working from model run
# bigc.fit <- fit.bigc
########## Load model data _____________________________________________________
bigc.modeldata <- readRDS(here("N_uptake_NEON/data/model_datalist/bigc.data.rds"))
########## Load full daily and hourly datasets _________________________________
# these include datatime info, and the model data does not)
# # 15-min dataset
# path <- here("N_uptake_NEON/data/neon_data_clean/bigc_clean.csv")
# bigc.df <- read_csv(path) %>%
# mutate(local_datetime = with_tz(local_datetime, tzone="US/Pacific"),
# model_datetime = local_datetime - hours(4))
# hourly dataset
path_h <- here("N_uptake_NEON/data/neon_data_clean/bigc_hourly_clean.csv")
bigc.df.h <- read_csv(path_h) %>%
mutate(local_datetime = with_tz(local_datetime, tzone="US/Pacific"),
model_datetime = local_datetime - hours(4)) First, check the fit for sigmas and betas (“sigma”, “sigma_U”, “b0”, “b1”), since these are notoriously tough to fit. Then put the spread and mean/ median/ quartile info for K, U, and N_e into a BIGC tibble (to be saved in one large csv at the end of the document). Finally, assemble tibbles for the visualizations, which include both modeled parameters and real data.
Inference for Stan model: anon_model.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.05 0.00 0.00 0.05 0.05 0.05 0.05 0.06 2811 1.00
sigma_U 0.16 0.00 0.04 0.09 0.13 0.16 0.18 0.23 129 1.02
b0 -9.15 0.03 0.74 -10.69 -9.62 -9.14 -8.67 -7.77 667 1.01
b1 0.82 0.00 0.08 0.66 0.77 0.82 0.88 1.00 683 1.01
Samples were drawn using NUTS(diag_e) at Tue Apr 1 09:45:41 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
########## N Model fit to data __________________________________________________
###### N model fit over time - predicted N vs measured N #######
N_Npred_ts_bigc <- plot_Nmodelfit_ts(data=N_output_bigc.df, start_jday = 35, end_jday = 60, plot_title = "NEON N data + model predictions: Upper Big Creek, CA")
# quartz(width=6,height=6)
N_Npred_ts_bigc### w. credible interval ribbon
N_Npred_tsCI_bigc <- plot_Nmodelfit_tsCI(data=N_output_bigc.df, start_jday = 35, end_jday = 60, plot_title = "NEON N data + model predictions: Upper Big Creek, CA")
# quartz(width=6,height=6)
N_Npred_tsCI_bigc###### predicted N vs measured N + 1:1 line #######
Npred_v_N_bigc <- plot_Npred_v_Nobs(data=N_output_bigc.df, plot_title="Predicted N vs. measured N, Upper Big Creek, CA")
# quartz(width=6, height=6)
Npred_v_N_bigc########## Equifinality checks _________________________________________________
####### Check equifinality!! K vs U #######
KvU_bigc <- plot_KvU(data=daily_output_bigc.df, plot_title = "Equifinality check, Upper Big Creek, CA: modeled K vs modeled U")
# quartz(width=7.5, height=6)
KvU_bigc####### Check equifinality!! K vs N_e #######
KvNe_bigc <- plot_KvNe(data=daily_output_bigc.df, plot_title = "Equifinality check, Upper Big Creek, CA: modeled K vs modeled equilibrium N")
# quartz(width=7.5, height=6)
KvNe_bigc####### Check equifinality!! U vs N_e #######
UvNe_bigc <- plot_UvNe(data=daily_output_bigc.df, plot_title = "Equifinality check, Upper Big Creek, CA: modeled U vs modeled equilibrium N")
# quartz(width=7.5, height=6)
UvNe_bigc########## Other parameter visualizations _______________________________________
####### K over time #######
K_v_time_bigc <- plot_K_ts(data=daily_output_bigc.df, plot_title = "modeled K over time, Upper Big Creek, CA")
# quartz(width=6.5, height=6)
K_v_time_bigc##### U over time #######
U_time_bigc <- plot_U_ts(data=daily_output_bigc.df, plot_title="Diel nitrate uptake (modeled), Upper Big Creek, CA")
# quartz(width=6, height=6)
U_time_bigc###### U vs sumlight #######
U_vs_light_bigc <- plot_UvLight(data=daily_output_bigc.df, plot_title = "Scatterplot of NO3 uptake and daily light: Upper Big Creek, CA")
# quartz(width=6.5, height=6)
U_vs_light_bigcWhat percentage of total uptake is the autotrophic uptake?
# NON-AUTOTROPHIC UPTAKE = Equilibrium nitrate * K
# TOTAL UPTAKE = U + (Equilibrium nitrate * K)
# Since U_auto is very small, there won't be much difference between 'non-autotrophic uptake' and 'total uptake'
daily_output_bigc.df <- daily_output_bigc.df %>%
mutate(U_het = N_e*K*z,
U_tot = U_het + U,
U_auto_perc = 100*U/U_tot)
U_auto_avg <- mean(daily_output_bigc.df$U_auto_perc) #11.7% on average
Uauto_perc_bigc <- plot_autoUperc(data=daily_output_bigc.df, plot_title = "Daily autotrophic uptake (as % of total uptake), Upper Big Creek, CA")
# quartz(width=7.5, height=6)
Uauto_perc_bigc #### Interactive w. ggplotly, but only gives the x,y info... maybe if I added datetime-associated labels?
# ggplotly(Uauto_plot_bigc)# path = here("N_uptake_NEON/data/model_output/bigc_output_daily")
# write_csv(daily_output_bigc.df, file=path)
#
# path.h = here("N_uptake_NEON/data/model_output/bigc_output_hourly")
# write_csv(N_output_bigc.df, file=path.h)Outputs and visualizations for Caribou Creek, Fairbanks North Star, AK
########## Load model fit ______________________________________________________
cari.fit <- readRDS(here("N_uptake_NEON/data/model_fits/cari.fit.rds"))
# if working from model run
# cari.fit <- fit.cari
########## Load model data _____________________________________________________
cari.modeldata <- readRDS(here("N_uptake_NEON/data/model_datalist/cari.data.rds"))
########## Load full hourly datasets __________________________________
# these include datatime info, and the model data does not)
# hourly dataset
path_h <- here("N_uptake_NEON/data/neon_data_clean/cari_hourly_clean.csv")
cari.df.h <- read_csv(path_h) %>%
mutate(local_datetime = with_tz(local_datetime, tzone="US/Alaska"),
model_datetime = local_datetime - hours(2))# 2-hr offset because it starts getting light at 02:45 in midsummer...First, check the fit for sigmas and betas (“sigma”, “sigma_U”, “b0”, “b1”), since these are notoriously tough to fit. Then put the spread and mean/ median/ quartile info for K, U, and N_e into a CARI tibble (to be saved in one large csv at the end of the document). Finally, assemble tibbles for the visualizations, which include both modeled parameters and real data.
Inference for Stan model: anon_model.
4 chains, each with iter=3000; warmup=1500; thin=1;
post-warmup draws per chain=1500, total post-warmup draws=6000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.10 0 0.00 0.10 0.10 0.10 0.11 0.11 10850 1
sigma_U 0.28 0 0.02 0.23 0.26 0.28 0.29 0.32 1579 1
b0 -2.26 0 0.29 -2.85 -2.45 -2.25 -2.06 -1.69 3517 1
b1 0.19 0 0.03 0.13 0.17 0.19 0.21 0.25 3516 1
Samples were drawn using NUTS(diag_e) at Tue Apr 15 15:00:35 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
########## N Model fit to data _________________________________________________
###### N model fit over time - predicted N vs measured N #######
N_Npred_ts_cari <- plot_Nmodelfit_ts(data=N_output_cari.df, start_jday = 145, end_jday = 155, plot_title = "NEON N data + model predictions: Caribou Creek, AK")
# quartz(width=6,height=6)
N_Npred_ts_cari### w. credible interval ribbon
N_Npred_tsCI_cari <- plot_Nmodelfit_tsCI(data=N_output_cari.df, start_jday = 145, end_jday = 155, plot_title = "NEON N data + model predictions: Caribou Creek, AK")
# quartz(width=6,height=6)
N_Npred_tsCI_cari###### predicted N vs measured N + 1:1 line #######
Npred_v_N_cari <- plot_Npred_v_Nobs(data=N_output_cari.df, plot_title="Predicted N vs. measured N, Caribou Creek, AK")
# quartz(width=6, height=6)
Npred_v_N_cari########## Equifinality checks ________________________________________________
####### Check equifinality!! K vs U #######
KvU_cari <- plot_KvU(data=daily_output_cari.df, plot_title = "Equifinality check, Caribou Creek, AK: modeled K vs modeled U")
# quartz(width=7.5, height=6)
KvU_cari####### Check equifinality!! K vs N_e #######
KvNe_cari <- plot_KvNe(data=daily_output_cari.df, plot_title = "Equifinality check, Caribou Creek, AK: modeled K vs modeled equilibrium N")
# quartz(width=7.5, height=6)
KvNe_cari####### Check equifinality!! U vs N_e #######
UvNe_cari <- plot_UvNe(data=daily_output_cari.df, plot_title = "Equifinality check, Caribou Creek, AK: modeled U vs modeled equilibrium N")
# quartz(width=7.5, height=6)
UvNe_cari########## Other parameter visualizations ______________________________________
####### K over time #######
K_v_time_cari <- plot_K_ts(data=daily_output_cari.df, plot_title = "modeled K over time, Caribou Creek, AK")
# quartz(width=6.5, height=6)
K_v_time_cari##### U over time #######
U_time_cari <- plot_U_ts(data=daily_output_cari.df, plot_title="Diel nitrate uptake (modeled), Caribou Creek, AK")
# quartz(width=6, height=6)
U_time_cari###### U vs sumlight #######
U_vs_light_cari <- plot_UvLight(data=daily_output_cari.df, plot_title = "Scatterplot of NO3 uptake and daily light: Caribou Creek, AK")
# quartz(width=6.5, height=6)
U_vs_light_cariWhat percentage of total uptake is the autotrophic uptake?
# NON-AUTOTROPHIC UPTAKE = Equilibrium nitrate * K
# TOTAL UPTAKE = U + (Equilibrium nitrate * K)
# Since U_auto is very small, there won't be much difference between 'non-autotrophic uptake' and 'total uptake'
daily_output_cari.df <- daily_output_cari.df %>%
mutate(U_het = N_e*K*z,
U_tot = U_het + U,
U_auto_perc = 100*U/U_tot)
U_auto_avg <- mean(daily_output_cari.df$U_auto_perc) #4.9% on average
Uauto_perc_cari <- plot_autoUperc(data=daily_output_cari.df, plot_title = "Daily autotrophic uptake (as % of total uptake), Caribou Creek, AK")
quartz(width=7.5, height=6)
Uauto_perc_cari
#### Interactive w. ggplotly, but only gives the x,y info... maybe if I added datetime-associated labels?
# ggplotly(Uauto_plot_cari)# path = here("N_uptake_NEON/data/model_output/cari_output_daily")
# write_csv(daily_output_cari.df, file=path)
#
# path.h = here("N_uptake_NEON/data/model_output/cari_output_hourly")
# write_csv(N_output_cari.df, file=path.h)Outputs and visualizations for Rio Cupeyes, San German Municipio, PR
########## Load model fit ______________________________________________________
cupe.fit <- readRDS(here("N_uptake_NEON/data/model_fits/cupe.fit.rds"))
# if working from model run
# cupe.fit <- fit.cupe
########## Load model data _____________________________________________________
cupe.modeldata <- readRDS(here("N_uptake_NEON/data/model_datalist/cupe.data.rds"))
########## Load full daily and hourly datasets __________________________________
# these include datatime info, and the model data does not)
# # 15-min dataset
# path <- here("N_uptake_NEON/data/neon_data_clean/cupe_clean.csv")
# cupe.df <- read_csv(path) %>%
# mutate(local_datetime = with_tz(local_datetime, tzone="America/Puerto_Rico"),
# model_datetime = local_datetime - hours(4))
# hourly dataset
path_h <- here("N_uptake_NEON/data/neon_data_clean/cupe_hourly_clean.csv")
cupe.df.h <- read_csv(path_h) %>%
mutate(local_datetime = with_tz(local_datetime, tzone="America/Puerto_Rico"),
model_datetime = local_datetime - hours(4)) First, check the fit for sigmas and betas (“sigma”, “sigma_U”, “b0”, “b1”), since these are notoriously tough to fit. Then put the spread and mean/ median/ quartile info for K, U, and N_e into a CUPE tibble (to be saved in one large csv at the end of the document). Finally, assemble tibbles for the visualizations, which include both modeled parameters and real data.
Inference for Stan model: anon_model.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.15 0.00 0.00 0.14 0.15 0.15 0.15 0.15 5813 1
sigma_U 0.24 0.00 0.01 0.22 0.24 0.24 0.25 0.26 2620 1
b0 -5.25 0.01 0.44 -6.08 -5.56 -5.25 -4.96 -4.37 2832 1
b1 0.56 0.00 0.05 0.46 0.53 0.56 0.60 0.66 2841 1
Samples were drawn using NUTS(diag_e) at Tue Apr 1 15:52:19 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
########## N Model fit to data __________________________________________________
###### N model fit over time - predicted N vs measured N #######
N_Npred_ts_cupe <- plot_Nmodelfit_ts(data=N_output_cupe.df, start_jday = 35, end_jday = 60, plot_title = "NEON N data + model predictions: Rio Cupeyes, PR")
# quartz(width=6,height=6)
N_Npred_ts_cupe### w. credible interval ribbon
N_Npred_tsCI_cupe <- plot_Nmodelfit_tsCI(data=N_output_cupe.df, start_jday = 35, end_jday = 60, plot_title = "NEON N data + model predictions: Rio Cupeyes, PR")
# quartz(width=6,height=6)
N_Npred_tsCI_cupe###### predicted N vs measured N + 1:1 line #######
Npred_v_N_cupe <- plot_Npred_v_Nobs(data=N_output_cupe.df, plot_title="Predicted N vs. measured N, Rio Cupeyes, PR")
# quartz(width=6, height=6)
Npred_v_N_cupe########## Equifinality checks _________________________________________________
####### Check equifinality!! K vs U #######
KvU_cupe <- plot_KvU(data=daily_output_cupe.df, plot_title = "Equifinality check, Rio Cupeyes, PR: modeled K vs modeled U")
# quartz(width=7.5, height=6)
KvU_cupe####### Check equifinality!! K vs N_e #######
KvNe_cupe <- plot_KvNe(data=daily_output_cupe.df, plot_title = "Equifinality check, Rio Cupeyes, PR: modeled K vs modeled equilibrium N")
# quartz(width=7.5, height=6)
KvNe_cupe####### Check equifinality!! U vs N_e #######
UvNe_cupe <- plot_UvNe(data=daily_output_cupe.df, plot_title = "Equifinality check, Rio Cupeyes, PR: modeled U vs modeled equilibrium N")
# quartz(width=7.5, height=6)
UvNe_cupe########## Other parameter visualizations ______________________________________
####### K over time #######
K_v_time_cupe <- plot_K_ts(data=daily_output_cupe.df, plot_title = "modeled K over time, Rio Cupeyes, PR")
# quartz(width=6.5, height=6)
K_v_time_cupe##### U over time #######
U_time_cupe <- plot_U_ts(data=daily_output_cupe.df, plot_title="Diel nitrate uptake (modeled), Rio Cupeyes, PR")
# quartz(width=6, height=6)
U_time_cupe###### U vs sumlight #######
U_vs_light_cupe <- plot_UvLight(data=daily_output_cupe.df, plot_title = "Scatterplot of NO3 uptake and daily light: Rio Cupeyes, PR")
# quartz(width=6.5, height=6)
U_vs_light_cupeWhat percentage of total uptake is the autotrophic uptake?
# NON-AUTOTROPHIC UPTAKE = Equilibrium nitrate * K * z
# TOTAL UPTAKE = U + (Equilibrium nitrate * K * z)
# Since U_auto is very small, there won't be much difference between 'non-autotrophic uptake' and 'total uptake'
daily_output_cupe.df <- daily_output_cupe.df %>%
mutate(U_het = N_e*K*z,
U_tot = U_het + U,
U_auto_perc = 100*U/U_tot)
U_auto_avg <- mean(daily_output_cupe.df$U_auto_perc) #4.7% on average
Uauto_perc_cupe <- plot_autoUperc(data=daily_output_cupe.df, plot_title = "Daily autotrophic uptake (as % of total uptake), Rio Cupeyes, PR")
# quartz(width=7.5, height=6)
Uauto_perc_cupe #### Interactive w. ggplotly, but only gives the x,y info... maybe if I added datetime-associated labels?
# ggplotly(Uauto_plot_cupe)# path = here("N_uptake_NEON/data/model_output/cupe_output_daily")
# write_csv(daily_output_cupe.df, file=path)
#
# path.h = here("N_uptake_NEON/data/model_output/cupe_output_hourly")
# write_csv(N_output_cupe.df, file=path.h)Outputs and visualizations for Pringle Creek, Wise County, TX
########## Load model fit ______________________________________________________
prin.fit <- readRDS(here("N_uptake_NEON/data/model_fits/prin.fit.rds"))
# if working from model run
# prin.fit <- fit.prin
########## Load model data _____________________________________________________
prin.modeldata <- readRDS(here("N_uptake_NEON/data/model_datalist/prin.data.rds"))
########## Load full daily and hourly datasets __________________________________
# these include datatime info, and the model data does not)
# # 15-min dataset
# path <- here("N_uptake_NEON/data/neon_data_clean/prin_clean.csv")
# prin.df <- read_csv(path) %>%
# mutate(local_datetime = with_tz(local_datetime, tzone="US/Central"),
# model_datetime = local_datetime - hours(4))
# hourly dataset
path_h <- here("N_uptake_NEON/data/neon_data_clean/prin_hourly_clean.csv")
prin.df.h <- read_csv(path_h) %>%
mutate(local_datetime = with_tz(local_datetime, tzone="US/Central"),
model_datetime = local_datetime - hours(4)) First, check the fit for sigmas and betas (“sigma”, “sigma_U”, “b0”, “b1”), since these are notoriously tough to fit. Then put the spread and mean/ median/ quartile info for K, U, and N_e into a PRIN tibble (to be saved in one large csv at the end of the document). Finally, assemble tibbles for the visualizations, which include both modeled parameters and real data.
Inference for Stan model: anon_model.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.31 0.00 0.00 0.30 0.30 0.31 0.31 0.32 6493 1
sigma_U 0.36 0.00 0.04 0.28 0.33 0.35 0.38 0.44 1104 1
b0 -6.48 0.03 1.26 -9.11 -7.28 -6.48 -5.65 -4.03 1608 1
b1 0.70 0.00 0.15 0.42 0.61 0.70 0.80 1.00 1621 1
Samples were drawn using NUTS(diag_e) at Wed Apr 16 00:03:30 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
########## N Model fit to data __________________________________________________
###### N model fit over time - predicted N vs measured N #######
N_Npred_ts_prin <- plot_Nmodelfit_ts(data=N_output_prin.df, start_jday = 90, end_jday = 120, plot_title = "NEON N data + model predictions: Pringle Creek, TX")
# quartz(width=6,height=6)
N_Npred_ts_prin### w. credible interval ribbon
N_Npred_tsCI_prin <- plot_Nmodelfit_tsCI(data=N_output_prin.df, start_jday = 90, end_jday = 120, plot_title = "NEON N data + model predictions: Pringle Creek, TX")
# quartz(width=6,height=6)
N_Npred_tsCI_prin###### predicted N vs measured N + 1:1 line #######
Npred_v_N_prin <- plot_Npred_v_Nobs(data=N_output_prin.df, plot_title="Predicted N vs. measured N, Pringle Creek, TX")
# quartz(width=6, height=6)
Npred_v_N_prin########## Equifinality checks _________________________________________________
####### Check equifinality!! K vs U #######
KvU_prin <- plot_KvU(data=daily_output_prin.df, plot_title = "Equifinality check, Pringle Creek, TX: modeled K vs modeled U")
# quartz(width=7.5, height=6)
KvU_prin####### Check equifinality!! K vs N_e #######
KvNe_prin <- plot_KvNe(data=daily_output_prin.df, plot_title = "Equifinality check, Pringle Creek, TX: modeled K vs modeled equilibrium N")
# quartz(width=7.5, height=6)
KvNe_prin####### Check equifinality!! U vs N_e #######
UvNe_prin <- plot_UvNe(data=daily_output_prin.df, plot_title = "Equifinality check, Pringle Creek, TX: modeled U vs modeled equilibrium N")
# quartz(width=7.5, height=6)
UvNe_prin########## Other parameter visualizations ______________________________________
####### K over time #######
K_v_time_prin <- plot_K_ts(data=daily_output_prin.df, plot_title = "modeled K over time, Pringle Creek, TX")
# quartz(width=6.5, height=6)
K_v_time_prin##### U over time #######
U_time_prin <- plot_U_ts(data=daily_output_prin.df, plot_title="Diel nitrate uptake (modeled), Pringle Creek, TX")
# quartz(width=6, height=6)
U_time_prin###### U vs sumlight #######
U_vs_light_prin <- plot_UvLight(data=daily_output_prin.df, plot_title = "Scatterplot of NO3 uptake and daily light: Pringle Creek, TX")
# quartz(width=6.5, height=6)
U_vs_light_prinWhat percentage of total uptake is the autotrophic uptake?
# NON-AUTOTROPHIC UPTAKE = Equilibrium nitrate * K * z
# TOTAL UPTAKE = U + (Equilibrium nitrate * K * z)
# Since U_auto is very small, there won't be much difference between 'non-autotrophic uptake' and 'total uptake'
daily_output_prin.df <- daily_output_prin.df %>%
mutate(U_het = N_e*K*z,
U_tot = U_het + U,
U_auto_perc = 100*U/U_tot)
U_auto_avg <- mean(daily_output_prin.df$U_auto_perc) #12.74% on average
Uauto_perc_prin <- plot_autoUperc(data=daily_output_prin.df, plot_title = "Daily autotrophic uptake (as % of total uptake), Pringle Creek, TX")
# quartz(width=7.5, height=6)
Uauto_perc_prin #### Interactive w. ggplotly, but only gives the x,y info... maybe if I added datetime-associated labels?
# ggplotly(Uauto_plot_prin)# path = here("N_uptake_NEON/data/model_output/prin_output_daily")
# write_csv(daily_output_prin.df, file=path)
#
# path.h = here("N_uptake_NEON/data/model_output/prin_output_hourly")
# write_csv(N_output_prin.df, file=path.h)Ouptuts and visualizations for West St Louis Creek, Grand, CO
########## Load model fit ______________________________________________________
wlou.fit <- readRDS(here("N_uptake_NEON/data/model_fits/wlou.fit.rds"))
# if working from model run
# wlou.fit <- fit.wlou
########## Load model data _____________________________________________________
# data from model
wlou.modeldata <- readRDS(here("N_uptake_NEON/data/model_datalist/wlou.data.rds"))
###### Load dataset (includes datatime info, which model data does not) _______
# # 15-min dataset
# path <- here("N_uptake_NEON/data/neon_data_clean/wlou_clean.csv")
# wlou.df <- read_csv(path) %>%
# mutate(local_datetime = with_tz(local_datetime, tzone="US/Mountain"),
# model_datetime = local_datetime - hours(4))
#
# hourly dataset
path_h <- here("N_uptake_NEON/data/neon_data_clean/wlou_hourly_clean.csv")
wlou.df.h <- read_csv(path_h) %>%
mutate(local_datetime = with_tz(local_datetime, tzone="US/Mountain"),
model_datetime = local_datetime - hours(4)) First, check the fit for sigmas and betas (“sigma”, “sigma_U”, “b0”, “b1”), since these are notoriously tough to fit. Then put the spread and mean/ median/ quartile info for K, U, and N_e into a WLOU tibble (to be saved in one large csv at the end of the document). Finally, assemble tibbles for the visualizations, which include both modeled parameters and real data.
Inference for Stan model: anon_model.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.05 0.00 0.00 0.05 0.05 0.05 0.05 0.05 5201 1
sigma_U 0.41 0.00 0.03 0.35 0.39 0.41 0.43 0.47 2253 1
b0 -2.64 0.01 0.57 -3.74 -3.03 -2.64 -2.25 -1.52 3282 1
b1 0.11 0.00 0.07 -0.02 0.07 0.11 0.16 0.24 3283 1
Samples were drawn using NUTS(diag_e) at Mon Mar 10 15:19:42 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
########## N Model fit to data ________________________________________________
###### N model fit over time #######
# N_and_Npred_wlou <- N_output_wlou.df %>%
# #group_by(year(model_datetime_wlou)) %>%
# filter(model_jday >= 35 & model_jday <= 60) %>% # to see the little boxes better...
# ggplot(aes(x=mod_hours)) +
# geom_point(aes(y=N_conc)) +
# geom_line(aes(y=N_conc_pred), col='red')+
# labs(
# x="Time (h)", y=expression("N"~(mmol~m^-3))
# ) +
# ggtitle("NEON: West St. Louis Creek N data + model predictions") +
# facet_wrap(~yr_jday) +
# theme_bw()
N_Npred_ts_wlou <- plot_Nmodelfit_ts(data=N_output_wlou.df, start_jday = 35, end_jday = 60, plot_title = "NEON N data + model predictions: West St. Louis Creek, CO ")
#quartz(width=6,height=6)
N_Npred_ts_wlou### w. credible interval ribbon
N_Npred_tsCI_wlou <- plot_Nmodelfit_tsCI(data=N_output_wlou.df, start_jday = 35, end_jday = 60, plot_title = "NEON N data + model predictions: West St. Louis Creek, CO")
# quartz(width=6,height=6)
N_Npred_tsCI_wlou###### N-pred vs N #######
# Npred_v_N_wlou <- ggplot(data = N_output_wlou.df, aes(x=N_conc, y=N_conc_pred)) +
# geom_point() +
# labs( x=expression("measured N"~(mu*mol~L^-1)), # fcn changed both to mmol~m^-3
# y=expression("modeled N"~(mu*mol~L^-1)) # ~ = a space, * = no space
# ) +
# ggtitle("Measured N vs predicted N, West St. Louis Creek, CO") +
# geom_abline(intercept = 0, slope = 1, color = "red") +
# theme_bw()
Npred_v_N_wlou <- plot_Npred_v_Nobs(data=N_output_wlou.df, plot_title="Predicted N vs. measured N, West St. Louis Creek, CO")
#quartz(width=6, height=6)
Npred_v_N_wlou########## Equifinality checks ________________________________________________
####### Check equifinality!! K vs U #######
# KvU_wlou <- ggplot(data=daily_output_wlou.df, aes(y=K_mod_avg_wlou, x=U_mod_avg_wlou)) +
# geom_point() +
# labs(y=expression("modeled K"~(d^-1)),
# x=expression("modeled U"~(mmol~m^-2~d^-1))
# ) +
# ggtitle("Equifinality check, West St. Louis Creek, CO: modeled K vs modeled U") +
# # geom_abline(intercept = 0, slope = 1, color = "red") + # doesn't show up - off the charts!
# theme_bw()
KvU_wlou <- plot_KvU(data=daily_output_wlou.df, plot_title = "Equifinality check, West St. Louis Creek, CO: modeled K vs modeled U")
#quartz(width=7.5, height=6)
KvU_wlou####### Check equifinality!! K vs N_e #######
# KvNe_wlou <- ggplot(data=daily_output_wlou.df, aes(y=K, x=N_e)) +
# geom_point() +
# labs(y=expression("modeled K"~(d^-1)),
# x=expression("modeled N_e"~(mmol~m^-3))
# ) +
# ggtitle("Equifinality check, West St. Louis Creek, CO: modeled K vs modeled equilibrium N") +
# # geom_abline(intercept = 0, slope = 1, color = "red") + # doesn't show up - off the charts!
# theme_bw()
KvNe_wlou <- plot_KvNe(data=daily_output_wlou.df, plot_title = "Equifinality check, West St. Louis Creek, CO: modeled K vs modeled equilibrium N")
#quartz(width=7.5, height=6)
KvNe_wlou####### Check equifinality!! U vs N_e #######
# UvNe_wlou <- ggplot(data=daily_output_wlou.df, aes(y=U, x=N_e)) +
# geom_point() +
# labs(y=expression("modeled U"~(mmol~m^-2~d^-1)),
# x=expression("modeled N_e"~(mmol~m^-3))
# ) +
# ggtitle("Equifinality check, West St. Louis Creek, CO: modeled U vs modeled equilibrium N") +
# theme_bw()
UvNe_wlou <- plot_UvNe(data=daily_output_wlou.df, plot_title = "Equifinality check, West St. Louis Creek, CO: modeled U vs modeled equilibrium N")
#quartz(width=7.5, height=6)
UvNe_wlou########## Other parameter visualizations ______________________________________
######## K over time #######
# K_v_time_wlou <- ggplot(data = daily_output_wlou.df, aes(x=d, y=K)) +
# geom_point() +
# # geom_point(y = sumlight.real, color = 'gold') +
# # ADD IN HIGH AND LOW CIs
# # xlab("Julian day") + ylab("modeled K (day -1)") + # daily change in N concentration
# labs( x= "Index",
# y=expression("modeled K"~(d^-1)) # ~ = a space, * = no space
# ) +
# #ylim(0,20) +
# ggtitle("modeled K over time, West St. Louis Creek, CO") +
# theme_bw()
K_v_time_wlou <- plot_K_ts(data=daily_output_wlou.df, plot_title = "modeled K over time, West St. Louis Creek, CO")
#quartz(width=6.5, height=6)
K_v_time_wlou##### to clip: no function yet
# K_time_clip_wlou <- ggplot(data = daily_output_wlou.df, aes(x=mod_day_wlou, y=K_mod_avg_wlou)) +
# geom_point() +
# # geom_point(y = sumlight.real, color = 'gold') +
# # ADD IN HIGH AND LOW CIs
# xlab("Julian day") + ylab("modeled K (day -1)") + # ??? UNITS?
# ylim(0,10) +
# ggtitle("modeled K over time, West St. Louis Creek, CO - ylim = 0,10") +
# theme_bw()
#
# quartz()
# K_time_clip_wlou
######## U over time #######
# U_time_wlou <- ggplot(data = daily_output_wlou.df, aes(x=d, y=U)) +
# geom_point() +
# # geom_point(y = sumlight.real_wlou, color = 'gold') +
# # ADD IN HIGH AND LOW CIs
# #xlab("Julian day") + ylab("modeled U (mmol/m2/day)") +
# labs( x= "Index",
# y=expression("modeled U"~(mmol~m^-2~d^-1)) # ~ = a space, * = no space
# ) +
# #ylim(0,1) +
# #ggtitle("modeled U over time, Big Creek 2019 pooled model w real light") +
# ggtitle("Diel nitrate uptake (modeled), West St. Louis Creek, CO") +
# theme_bw()
U_time_wlou <- plot_U_ts(data=daily_output_wlou.df, plot_title="Diel nitrate uptake (modeled), West St. Louis Creek, CO")
#quartz(width=6, height=6)
U_time_wlou# U_time_clip_wlou <- ggplot(data = daily_output_wlou.df, aes(x=mod_day_wlou, y=U_mod_avg_wlou)) +
# geom_point() +
# # geom_point(y = sumlight.real, color = 'gold') +
# # ADD IN HIGH AND LOW CIs
# xlab("Julian day") + ylab("modeled U (mmol/m2/day)") +
# ylim(0,1) +
# ggtitle("modeled U over time, West St. Louis Creek, CO - ylim = 0-1") +
# theme_bw()
#
# quartz()
# U_time_clip_wlou
######## U vs sumlight #######
# U_vs_light_wlou <- ggplot(data = daily_output_wlou.df, aes(x=sumlightReal, y=U)) +
# geom_point() +
# #xlab("true light (satellite)") + ylab("modeled U (mmol/m2/day)") +
# # xlab("light (satellite)") + ylab("modeled U (mmol/m2/day)") +
# labs( x= "light (satellite)",
# y=expression("modeled U"~(mmol~m^-2~d^-1)) # ~ = a space, * = no space
# ) +
# scale_x_log10() + scale_y_log10() +
# ggtitle("Scatterplot of NO3 uptake and daily light: West St. Louis Creek, CO") +
# #ylim = c(-0.2, 1) +
# theme_bw()
U_vs_light_wlou <- plot_UvLight(data=daily_output_wlou.df, plot_title = "Scatterplot of NO3 uptake and daily light: West St. Louis Creek, CO")
# quartz(width=6.5, height=6)
U_vs_light_wlou# NON-AUTOTROPHIC UPTAKE = Equilibrium nitrate * K
# TOTAL UPTAKE = U + (Equilibrium nitrate * K)
# Since U_auto is very small, there won't be much difference between 'non-autotrophic uptake' and 'total uptake'
daily_output_wlou.df <- daily_output_wlou.df %>%
mutate(U_het = N_e*K*z,
U_tot = U_het + U,
U_auto_perc = 100*U/U_tot)
U_auto_avg <- mean(daily_output_wlou.df$U_auto_perc) #10.3% on average
# Uauto_plot_wlou <- ggplot(data=daily_output_wlou.df, aes(x=d, y=U_auto_perc)) +
# geom_point() +
# xlab("Index") + ylab("Percent autotrophic uptake") +
# ggtitle("Percent of autotrophic uptake each day, West St. Louis Creek, CO") +
# theme_bw()
Uauto_perc_wlou <- plot_autoUperc(data=daily_output_wlou.df, plot_title = "Daily autotrophic uptake (as % of total uptake), West St. Louis Creek, CO")
# quartz(width=7.5, height=6)
Uauto_perc_wlou #### Interactive w. ggplotly, but only gives the x,y info... maybe if I added datetime-associated labels?
# ggplotly(Uauto_plot_wlou)# path = here("N_uptake_NEON/data/model_output/wlou_output_daily")
# write_csv(daily_output_wlou.df, file=path)
#
# path.h = here("N_uptake_NEON/data/model_output/wlou_output_hourly")
# write_csv(N_output_wlou.df, file=path.h)Combine the site tibbles into one large tibble, then save it
########## Reload all datasets as needed _______________________________________
########## Bind then save daily parameters (K, U, N_e) _________________________
# daily_summary_all.df <- bind_rows()
# path <- here()
# saveRDS(daily_summary_all.df, path)
########## Bind then save hourly parameters (N_pred) ___________________________
# hourly_summary_all.df <- bind_rows()
# path2 <- here()
# saveRDS(hourly_summary_all.df, path2)########## Violin plots _______________________________________________
##### K
##### U
##### N_e